home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_6.3 / XSGT / XSGT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  28.9 KB  |  1,101 lines
  1. /* 
  2.  * xsgt.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * (e)X(tract) S(ite) G(eometry and) T(opology) 
  11.  *
  12.  * routine to merge geometry and topology data for bubble patterns from:
  13.  * -->files of type .gdt, generated by xah
  14.  * -->files of type .std, generated by xvora 
  15.  * output file: .sgt 
  16.  *
  17.  * various geometrical and topological quantities are evaluated:
  18.  * -->p(n, A), joint coordination number, area distribution function
  19.  * -->s_ij, distance of closest approach
  20.  *
  21.  */
  22.  
  23. #include "sgt.h"
  24.  
  25. #define MAX_RECORD_SIZE    600
  26.  
  27.  
  28. #define LINE_LEN    96          /* max length of line in data file */
  29. #define FN_LEN      16          /* length of data file name */
  30.  
  31. #define ON      1
  32. #define OFF         0
  33. #define ECHO_FILE_NAME  ON
  34.  
  35. #define GDT_DEBUG
  36. #undef  SHOW_HSORT
  37. #undef  FRMS_DEBUG
  38. #undef  KNNS_ARRAY_VER          /* see comments in knns.c */
  39.  
  40.  
  41. /*
  42.  * global variables
  43.  */
  44. extern char *optarg;
  45. extern int optind, opterr;
  46.  
  47. int READ_AFLAG = 0;
  48. int READ_GFLAG = 0;
  49. int CHECK_EDGE = 1;
  50. int WRITE_FILE = 0;
  51.  
  52.  
  53. #define BS_HIST     "\nhistogram for boundary-separations"
  54. #define BIN_NUMBER  51          /* no bins in histogram */
  55. #define NSIJ_MIN    0.0         /* min, max for norm bound sep */
  56. #define NSIJ_MAX    2.5
  57.  
  58. char *HIST_HEADER;
  59. int HIST_TYPE = -1;
  60. int N_BINS = BIN_NUMBER;
  61. int MIN_SET = 0;
  62. int MAX_SET = 0;
  63.  
  64. int EVAL_KNN = 0;
  65. int SHOW_INPUT = 0;
  66. int SHOW_ALL = 1;
  67. int CONTRACT = 0;
  68.  
  69.  
  70. int KNN_MAX = 2;                /* max nof NN shells to be constructed */
  71.  
  72.  
  73. /*
  74.  * compare coordinates of two sites: allow for deviations of +/- 1;
  75.  */
  76. int
  77. cmp_site_coord (float x1, float y1, float x2, float y2)
  78. {
  79.  
  80.   if (F_TO_INT (fabs (x1 - x2)) <= 1) {
  81.     if (F_TO_INT (fabs (y1 - y2)) <= 1)
  82.       return (1);
  83.   }
  84.   return (-1);
  85. }
  86.  
  87.  
  88. /*
  89.  * comparison function for qsort():
  90.  * sort array of struct kNNShell in order of increasing values
  91.  * of the area of the root Site, that is, the value of the 0-th
  92.  * element in the array aka[];
  93.  */
  94. int
  95. cmp_root_a (t1, t2)
  96. /*const void *t1, *t2; */
  97.      struct kNNShell *t1, *t2;
  98. {
  99.   struct kNNShell *sk1, *sk2;
  100.  
  101.   sk1 = ((struct kNNShell *) t1);
  102.   sk2 = ((struct kNNShell *) t2);
  103.  
  104.   return ((int) SIGN (sk1->aka[0] - sk2->aka[0]));
  105. }
  106.  
  107.  
  108. /*
  109.  * comparison function for qsort():
  110.  * sort array of tuples in order of increassing x-values
  111.  */
  112. int
  113. compare (t1, t2)
  114. /*const void *t1, *t2; */
  115.      float *t1, *t2;
  116. {
  117.   float f1, f2;
  118.  
  119.   f1 = *((float *) t1);
  120.   f2 = *((float *) t2);
  121.   return ((int) SIGN (f1 - f2));
  122. }
  123.  
  124.  
  125. /*
  126.  * vprintf() functions:
  127.  * write to std output and, if option WRITE_FILE set, to file wbuf (stream)
  128.  */
  129. void
  130. gprintf (FILE * fpOut, char *fmt,...)
  131. {
  132.   va_list arg_ptr;
  133.  
  134.   va_start (arg_ptr, fmt);
  135.   vprintf (fmt, arg_ptr);
  136.   if (WRITE_FILE == 1)
  137.     vfprintf (fpOut, fmt, arg_ptr);
  138.   va_end (arg_ptr);
  139.  
  140. }
  141.  
  142.  
  143.  
  144. void
  145. fail_alloc (char *str, int code)
  146. {
  147.   printf ("\n...memory alloc for %s failed\n", str);
  148.   exit (code);
  149. }
  150.  
  151.  
  152. /*
  153.  * error message
  154.  */
  155. void
  156. exitmess (char *prompt, int status)
  157. {
  158.   printf (prompt);
  159.   printf ("\n");
  160.  
  161.   exit (status);
  162. }
  163.  
  164.  
  165. void
  166. usage (char *progname)
  167. {
  168.   progname = last_bs (progname);
  169.   printf ("USAGE: %s infile [-k S] [-n n] [-e] [-w] [-L]\n", progname);
  170.   printf ("\n%s merges geometry and topology data for droplet patterns\n", progname);
  171.   printf ("based on previous analysis by xah and xvora\n\n");
  172.   printf ("ARGUMENTS:\n");
  173.   printf ("        infile: input filename of type .xxx (current default)\n");
  174.   printf ("                NOTE: full names of the two required input files\n");
  175.   printf ("                infile.gdt and infile.std, and optionally an\n");
  176.   printf ("                output file, are generated after stripping the\n");
  177.   printf ("                .xxx suffix; input data file types:\n");
  178.   printf ("                     .gdt, generated by xah\n");
  179.   printf ("                     .std, generated by xvora\n\n");
  180.   printf ("OPTIONS:\n");
  181.   printf ("   -k S: construct k-NN shells, k<=S (default=%d)\n", KNN_MAX);
  182.   printf ("   -n n: set number of bins to n (default=%d)\n", N_BINS);
  183.   printf ("     -w: write output file of type .sgt (site geom topol)\n");
  184.   printf ("     -L: print Software License for this module\n");
  185.   exit (1);
  186. }
  187.  
  188.  
  189. /*
  190.  * skip n lines
  191.  */
  192. void
  193. skip_n_lines (file_pointer, n_lines_to_skip)
  194.      FILE *file_pointer;
  195.      int n_lines_to_skip;
  196. {
  197.   int i;
  198.   char ch;
  199.  
  200.   for (i = 0; i < n_lines_to_skip; i++) {
  201.     while ((ch = getc (file_pointer)) != '\n');
  202.   }
  203. }
  204.  
  205.  
  206. /*
  207.  * contract original Site array, by checking out-of-bounds flags:
  208.  * eliminate all Sites which have been flagged
  209.  */
  210. int
  211. contract_vsa (struct vSite **cvsa, struct vSite *vsa, int ns)
  212. {
  213.   int is, nas;
  214.   struct vSite *vS;
  215.   nas = 0;
  216.   for (is = 0; is < ns; is++) {
  217.     if (((vS = (vsa + is))->eFlag == UnSet) &&
  218.         (vS->aoiFlag == UnSet)) {
  219.       cvsa[is] = (vsa + is);
  220.       nas++;
  221.     }
  222.   }
  223.   return (nas);
  224. }
  225.  
  226.  
  227. void
  228. main (int argc, char **argv)
  229. {
  230.   int i_arg;
  231.  
  232.   int i, ich, is, k;
  233.   int inn, nnn, nnn_max, ns;
  234.   int sitenbr;
  235.   int n_tdt;
  236.   int eFlag, a_flag;
  237.  
  238.   char colon, komma, lpar, rpar;
  239.  
  240.   static char *inp_suffix =
  241.   {".xxx"};                     /* dummy inp file suffix */
  242.   static char *inp_asuffix =
  243.   {".gdt"};                     /* default inp file suffix */
  244.   static char *inp_gsuffix =
  245.   {".std"};                     /* default inp file suffix */
  246.   static char *xxx_type =
  247.   {".xxx"};
  248.   static char *wsuffix =
  249.   {".sgt"};                     /* suffix for output file name */
  250.   char *rbuf;
  251.  
  252.   static char rabuf[32], rgbuf[32];
  253.   static char wbuf[20];
  254.   char buf[32], ln_buf[LINE_LEN], fn_buf[16];
  255.   int stop_file_io = -1;
  256.  
  257.   int ngdt, np, nep;
  258.   float *a, *xs, *ys, *area, *sig;
  259.   float xc, yc;
  260.   struct vSite *vsa, **cvsa;
  261.  
  262.   FILE *fpaIn, *fpgIn, *fpOut;
  263.  
  264.   int idn;
  265.   int *ondn, *ndn;              /* arrays cont nof domains with n NN */
  266.   int nbs;
  267.   float **bs;                   /* ar of pntrs to bdy sep betw NN pairs */
  268.   unsigned int **a_n;           /* ar of pntrs to area of n-fold coord bubs */
  269.   double *a_nbar;               /* mean area for bubs of coord no n */
  270.   double del_an, *sig_an;       /* stand dev corresp to a_nbar */
  271.   double ma_tdt, mr_tdt, ma_gdt, mva;
  272.   struct Pix *dpa;
  273.  
  274.   int ih, nh;
  275.   float bw, min = (float) -1.0, max = (float) -1.0;
  276.   double mean, sdev;
  277.   float *hbs, *ha;              /* histogram of bs values */
  278.  
  279.  
  280.   float **m_n;                  /* ar of pntrs to avrg nof NN of NN of n-fold site */
  281.   double *m_nbar;               /* mean values of m_n, as func of n */
  282.   double del_mn, *sig_mn;       /* stand dev corresp to m_nbar */
  283.   int nas;
  284.  
  285.  
  286.   /* list structs */
  287.   struct kNNShell *sk;          /* array of pntrs to lists of kNN shells */
  288.   struct linklist *sha;         /* array of kNN shells (lists) */
  289.  
  290.   double frms;
  291.  
  292. /* 
  293.  * cmd line options:
  294.  */
  295.   static char *optstring = "r:k:n:ew";
  296.  
  297.   if (argc < 2)
  298.     usage (argv[0]);
  299.  
  300.  
  301.   rbuf = argv[1];
  302.  
  303.   /* construct input file names */
  304.   ich = 0;
  305.   while ((*(rabuf + ich) = *(rgbuf + ich) = *(rbuf + ich)) != '.')
  306.     ich++;
  307.   for (is = 0; is < 4; is++) {
  308.     *(rabuf + (ich) + is) = *(inp_asuffix + is);
  309.     *(rgbuf + (ich) + is) = *(inp_gsuffix + is);
  310.   }
  311.  
  312.   if ((fpaIn = fopen (rabuf, "r")) == NULL) {
  313.     printf ("\n...cannot open input file %s\n", rabuf);
  314.     exit (1);
  315.   }
  316.   nep = 3;
  317.   READ_AFLAG = 1;
  318.  
  319.   if ((fpgIn = fopen (rgbuf, "r")) == NULL) {
  320.     printf ("\n...cannot open input file %s\n", rgbuf);
  321.     exit (1);
  322.   }
  323.   READ_GFLAG = 1;
  324.  
  325.  
  326.   printf ("read data from files %s and %s\n", rabuf, rgbuf);
  327.   /* construct output file name */
  328.   ich = 0;
  329.   while ((*(wbuf + ich) = *(rbuf + ich)) != '.')
  330.     ich++;
  331.   for (is = 0; is < 4; is++)
  332.     *(wbuf + (ich) + is) = *(wsuffix + is);
  333.  
  334.   /* strip input file suffix */
  335.   for (is = 0; is < 4; is++)
  336.     *(inp_suffix + is) = *(rbuf + (ich) + is);
  337.  
  338.   if (strcmp (inp_suffix, xxx_type) != 0) {
  339.     printf ("\n...unknown input file type encountered:");
  340.     printf (" default: %s\n", xxx_type);
  341.     fclose (fpgIn);
  342.     fclose (fpaIn);
  343.     exit (1);
  344.   }
  345.  
  346. /*
  347.  * parse command line
  348.  */
  349.   optind = 2;
  350.   opterr = ON;                  /* give error messages */
  351.  
  352.  
  353.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  354.     switch (i_arg) {
  355.     case 'k':
  356.       printf ("...option %c: ", i_arg);
  357.       printf ("construct %d-NN shells for all Sites\n",
  358.               KNN_MAX = atoi (optarg));
  359.       EVAL_KNN = ON;
  360.       break;
  361.  
  362.     case 'n':
  363.       printf ("...option %c: ", i_arg);
  364.       printf ("set number of bins to %d\n",
  365.               N_BINS = atoi (optarg));
  366.       break;
  367.  
  368.     case 'e':
  369.       printf ("...option %c: ", i_arg);
  370.       printf ("disable Site status (eFlag, aoiFlag) check\n");
  371.       CHECK_EDGE = OFF;
  372.       break;
  373.  
  374.     case 'w':
  375.       printf ("\n...option %c: ", i_arg);
  376.       if ((fpOut = fopen (wbuf, "w")) == NULL) {
  377.         printf ("\n...cannot open %s for writing\n",
  378.                 wbuf);
  379.         exit (1);
  380.       }
  381.       printf ("write output to file %s\n", wbuf);
  382.       WRITE_FILE = 1;
  383.       break;
  384.  
  385.     case 'L':
  386.       print_sos_lic ();
  387.       exit (0);
  388.  
  389.     default:
  390.       printf ("\n...unknown cmd line argument\n");
  391.       exit (1);
  392.       break;
  393.     }
  394.   }
  395.   if ((READ_AFLAG == 0) || (READ_GFLAG == 0)) {
  396.     printf ("\n...must invoke option -r to read input files!!\n");
  397.     exit (1);
  398.   }
  399.  
  400. /* 
  401.  * merge .gdt and .std data files 
  402.  */
  403.  
  404. /*
  405.  * retrieve data from .gdt file (same type as 2-column standard file)
  406.  */
  407.   if ((ngdt = get_record_size (fpaIn)) > MAX_RECORD_SIZE) {
  408.     printf ("\n...record size of %d exceeds limit of %d\n",
  409.             ngdt, MAX_RECORD_SIZE);
  410.     exit (1);
  411.   }
  412.   if ((np = get_prm_size (fpaIn)) != nep) {
  413.     printf ("\n...n_parms = %d! expect: n_parms = %d\n",
  414.             np, nep);
  415.     exit (1);
  416.   }
  417.  
  418. /*
  419.  * mem alloc
  420.  */
  421.   if ((a = (float *) calloc (np, sizeof (float))) == NULL)
  422.       fail_alloc ("a", 1);
  423.  
  424.   if ((xs = (float *) calloc (ngdt, sizeof (float))) == NULL)
  425.       fail_alloc ("xs", 1);
  426.   if ((ys = (float *) calloc (ngdt, sizeof (float))) == NULL)
  427.       fail_alloc ("ys", 1);
  428.   if ((area = (float *) calloc (ngdt, sizeof (float))) == NULL)
  429.       fail_alloc ("area", 1);
  430.   if ((sig = (float *) calloc (ngdt, sizeof (float))) == NULL)
  431.       fail_alloc ("sig", 1);
  432.  
  433.  
  434. /* fetch data */
  435.   get_col3_data (fpaIn, np, a, ngdt, xs, ys, area, sig);
  436.   fclose (fpaIn);
  437.  
  438. #ifdef GDT_DEBUG
  439.   printf ("\n...data from input file %s:\n", rabuf);
  440.   printf ("\n...ngdt: %d  n_parms: %d\n", ngdt, np);
  441.  
  442.   for (i = 0; i < np; i++)
  443.     printf ("a[%d] = %f", i, *(a + i));
  444.   printf ("\n");
  445. #endif
  446.  
  447.   ma_gdt = 0.0;
  448.   for (i = 0; i < ngdt; i++) {
  449. #ifdef GDT_DEBUG
  450.     printf ("xs[%d] = %f, ys[%d] = %f, area[%d] = %f\n",
  451.             i, *(xs + i), i, *(ys + i), i, *(area + i));
  452. #endif
  453.     ma_gdt += (double) (*(area + i));
  454.   }
  455.   ma_gdt /= (double) ngdt;
  456.   printf ("\n...mean area ( .gdt): %lf\n", ma_gdt);
  457.  
  458.  
  459.   if ((vsa = (struct vSite *) calloc (ngdt, sizeof (struct vSite))) == NULL)
  460.       fail_alloc ("vsa", 1);
  461.  
  462. /*
  463.  * retrieve data from .std file
  464.  */
  465.   skip_n_lines (fpgIn, 44);
  466.   fgets (ln_buf, LINE_LEN, fpgIn);
  467.   sscanf (ln_buf, "%[^:]%c%s", buf, &colon, fn_buf);
  468.   printf ("\n...data from input file %s\n", fn_buf);
  469.  
  470.   fgets (ln_buf, LINE_LEN, fpgIn);
  471.   sscanf (ln_buf, "%[^:]%c%d", buf, &colon, &n_tdt);
  472.   printf ("...total nof site entries: %d\n", n_tdt);
  473.  
  474. /*
  475.  * the following code handles .std files containing just the site
  476.  * information (followed by EOF) as well as those with trailing
  477.  * histogram entries
  478.  */
  479. /*
  480.  * first pass through .std file:
  481.  * -->determine number of sites: should be identical to n_rec
  482.  * -->determine nnn_max;
  483.  */
  484.   ns = 0;
  485.   sitenbr = nnn = nnn_max = 0;
  486.   stop_file_io = -1;
  487.  
  488.   skip_n_lines (fpgIn, 2);
  489.   while (stop_file_io != 1) {
  490.     if (fgets (ln_buf, LINE_LEN, fpgIn) != (char *) NULL) {
  491.       sscanf (ln_buf,
  492.               "%[^:]%c%d%[^(]%c%f%c%f%[^:]%c%d%[^:]%c%d%[^:]%c%d",
  493.               buf, &colon, &sitenbr,
  494.               buf, &lpar, &xc, &komma, &yc,
  495.               buf, &colon, &eFlag,
  496.               buf, &colon, &a_flag,
  497.               buf, &colon, &nnn);
  498.  
  499.       if (sitenbr == n_tdt - 1)
  500.         stop_file_io = 1;
  501.  
  502.       if (nnn > nnn_max)
  503.         nnn_max = nnn;
  504.       skip_n_lines (fpgIn, 1 + 1 + nnn + 2);
  505.       ns++;
  506.     }
  507.     else {
  508.       exitmess ("ERROR reading .std file\n", 1);
  509.     }
  510.   }
  511.   fclose (fpgIn);
  512.   printf ("...first pass through data found:");
  513.   printf ("   %d sites, max nnn: %d\n", ns, nnn_max);
  514.  
  515.   if (ns != ngdt) {
  516.     printf ("WARNING: ");
  517.     printf (" sz of .gdt file, %d, does not match sz of .std file, %d\n",
  518.             ngdt, ns);
  519.  
  520.     exitmess ("-->remedy discrepancy before proceeding\n", 1);
  521.  
  522.   }
  523.  
  524.   for (is = 0; is < ns; is++) {
  525.     if (((vsa + is)->nnd = (float *) calloc (nnn_max, sizeof (float))) == NULL)
  526.         fail_alloc ("(vsa+is)->nnd", 1);
  527.  
  528.  
  529.     if (((vsa + is)->nnsi = (int *) calloc (nnn_max, sizeof (int))) == NULL)
  530.         fail_alloc ("(vsa+is)->nnsi", 1);
  531.  
  532. /* initialize arrays within structures */
  533.     for (inn = 0; inn < nnn_max; inn++) {
  534.       (vsa + is)->nnd[inn] = (float) -1.0;
  535.       (vsa + is)->nnsi[inn] = (int) -1;
  536.     }
  537.  
  538.     if ((ondn = (int *) calloc (nnn_max + 1, sizeof (int))) == NULL)
  539.         fail_alloc ("ondn", 1);
  540.     if ((ndn = (int *) calloc (nnn_max + 1, sizeof (int))) == NULL)
  541.         fail_alloc ("ndn", 1);
  542.  
  543.   }
  544.  
  545. /*
  546.  * second pass: 
  547.  * -->scan .std input file, 
  548.  * -->match site coords with those found in .gdt file
  549.  */
  550.   if ((fpgIn = fopen (rgbuf, "r")) == NULL)
  551.     exitmess ("cannot open input file", 1);
  552.   skip_n_lines (fpgIn, 44);
  553.  
  554.   skip_n_lines (fpgIn, 2 + 2);
  555.   ns = 0;
  556.   stop_file_io = -1;
  557.   while (stop_file_io != 1) {
  558.     if (fgets (ln_buf, LINE_LEN, fpgIn) != (char *) NULL) {
  559.       sscanf (ln_buf,
  560.               "%[^:]%c%d%[^(]%c%f%c%f%[^:]%c%d%[^:]%c%d%[^:]%c%d",
  561.               buf, &colon, &((vsa + ns)->sitenbr),
  562.               buf, &lpar, &((vsa + ns)->coord.x),
  563.               &komma, &((vsa + ns)->coord.y),
  564.               buf, &colon, &((vsa + ns)->eFlag),
  565.               buf, &colon, &((vsa + ns)->aoiFlag),
  566.               buf, &colon, &((vsa + ns)->nnn));
  567.  
  568.       if ((sitenbr = (vsa + ns)->sitenbr) == n_tdt - 1)
  569.         stop_file_io = 1;
  570.  
  571.       ondn[(vsa + ns)->nnn]++;
  572.  
  573.       fgets (ln_buf, LINE_LEN, fpgIn);
  574.       sscanf (ln_buf, "%[^:]%c%u",
  575.               buf, &colon, &((vsa + ns)->v_area));
  576.  
  577.       skip_n_lines (fpgIn, 1);
  578.       for (inn = 0; inn < (vsa + ns)->nnn; inn++) {
  579.         fgets (ln_buf, LINE_LEN, fpgIn);
  580.         sscanf (ln_buf, "%[^(]%c%d%c%f%c",
  581.                 buf, &lpar, &((vsa + ns)->nnsi[inn]),
  582.                 &colon, &((vsa + ns)->nnd[inn]), &rpar);
  583.       }
  584.       skip_n_lines (fpgIn, 2);
  585.  
  586.  
  587. /*
  588.  * match site coords and extract area value (assumes y- and x-ordered data)
  589.  */
  590.       if (cmp_site_coord (*(xs + ns), *(ys + ns),
  591.                           (float) (vsa + ns)->coord.x,
  592.                           (float) (vsa + ns)->coord.y) == 1) {
  593.         (vsa + ns)->area = (unsigned int) (*(area + ns));
  594.       }
  595.       else {
  596.         printf ("\n...no match: ");
  597.         printf ("(xs=%5.2f, ys=%5.2f) ",
  598.                 *(xs + ns), *(ys + ns));
  599.         printf ("(vsx=%5.2f, vsy=%5.2f)\n",
  600.                 (vsa + ns)->coord.x, (vsa + ns)->coord.y);
  601.         exitmess ("-->something wrong", 1);
  602.       }
  603.       ns++;
  604.     }
  605.   }
  606.   fclose (fpgIn);
  607.  
  608.  
  609. /*
  610.  * write output
  611.  */
  612.   if (WRITE_FILE == 1) {
  613.     gprintf (fpOut, "%d %d\n", ngdt, nnn_max);
  614.     gprintf (fpOut, "nof domains with coord no n=2(1)nnn_max: ");
  615.     for (inn = 2; inn <= nnn_max; inn++) {
  616.       gprintf (fpOut, " %d ", *(ondn + inn));
  617.     }
  618.     gprintf (fpOut, "\n");
  619.  
  620.     for (is = 0; is < ns; is++) {
  621.       gprintf (fpOut, "%3d %4.1f %4.1f %u %u %d %d %2d\n",
  622.                (vsa + is)->sitenbr,
  623.                (vsa + is)->coord.x, (vsa + is)->coord.y,
  624.                (vsa + is)->v_area, (vsa + is)->area,
  625.                (vsa + is)->eFlag, (vsa + is)->aoiFlag,
  626.                (vsa + is)->nnn);
  627.  
  628.       for (inn = 0; inn < (vsa + is)->nnn; inn++) {
  629.         gprintf (fpOut, "  %5d  %5.1f\n",
  630.                  (vsa + is)->nnsi[inn], (vsa + is)->nnd[inn]);
  631.       }
  632.     }
  633.   }
  634.  
  635.  
  636. /*
  637.  * contract Site array vsa by eliminating all Sites which have been 
  638.  * flagged out-of-bounds (eFlag, aoiFlag)
  639.  *
  640.  * following code (CONTRACT == 1) is largely untested: it might be useful
  641.  * to perform this contraction to avoid having to check, as done currently 
  642.  * for all later functions (in anal_gt.c, knns.c, etc) whether original Sites
  643.  * have been flagged out-of-bounds; this checking eliminates such Sites as
  644.  * reference (or ``root'') Sites, but leaves them present as possible kNN
  645.  * of other reference Sites, while contraction would remove them altogether
  646.  *
  647.  * note: performing the Site checking affects the counts of nnn, maintained
  648.  *   in the index array ndn: entries will generally differ from those 
  649.  *   in ondn; to make subsequent array alloc more robust (independent
  650.  *   of whether or not Site checking was enforced in pertinent functions
  651.  *   (e.g. p_of_a(), m_n(), etc), arrays are slightly over-dimensioned
  652.  *   by relying on ondn, the original entries in ndn with all Sites
  653.  *   admitted;
  654.  */
  655.   if (CONTRACT == 1) {
  656.     if ((cvsa = (struct vSite **) calloc (ns, sizeof (struct vSite *))) == NULL)
  657.         fail_alloc ("cvsa", 1);
  658.  
  659.     ns = contract_vsa (cvsa, vsa, ns);
  660.     printf ("\n...Site array contraction leaves %d active Sites\n", ns);
  661.  
  662.     if ((cvsa = (struct vSite **) realloc (cvsa, ns * sizeof (struct vSite *))) ==
  663.         NULL)
  664.         fail_alloc ("cvsa(realloc)", 1);
  665.  
  666.  
  667.     for (inn = 0; inn < nnn_max; inn++)
  668.       ndn[inn] = 0;
  669.     nnn_max = 0;
  670.     for (is = 0; is < ns; is++) {
  671.       if (cvsa[is]->nnn > nnn_max)
  672.         nnn_max = cvsa[is]->nnn;
  673.       ndn[cvsa[is]->nnn]++;
  674.     }
  675.     printf ("%d %d\n", ns, nnn_max);
  676.     printf ("nof domains with coord no n=2(1)nnn_max: ");
  677.     for (inn = 2; inn <= nnn_max; inn++) {
  678.       printf (" %d ", *(ndn + inn));
  679.     }
  680.     printf ("\n");
  681.  
  682.     for (is = 0; is < ns; is++) {
  683.       printf ("%3d %4.1f %4.1f %u %2d\n",
  684.               cvsa[is]->sitenbr,
  685.               cvsa[is]->coord.x, cvsa[is]->coord.y,
  686.               cvsa[is]->area, cvsa[is]->nnn);
  687.  
  688.       for (inn = 0; inn < cvsa[is]->nnn; inn++) {
  689.         printf ("  %5d  %5.1f\n",
  690.                 cvsa[is]->nnsi[inn], cvsa[is]->nnd[inn]);
  691.       }
  692.     }
  693.   }
  694.  
  695.  
  696. /*
  697.  * extract various geometrical and topological quantities 
  698.  * describing bubble domain patterns
  699.  */
  700.  
  701. /* 
  702.  * joint probability distribution of coordination number, n, and area, A 
  703.  *
  704.  * --> examine number of NNs as a function of domain area (->Lewis' law)
  705.  */
  706.  
  707.   if ((a_n = (unsigned int **) calloc (nnn_max, sizeof (unsigned int *))) == NULL)
  708.       fail_alloc ("a_n", 1);
  709.  
  710.   for (inn = 0; inn < nnn_max; inn++) {
  711.     if ((a_n[inn] = (unsigned int *) calloc (ondn[inn] + 1, sizeof (unsigned
  712.                                                              int))) == NULL)
  713.         fail_alloc ("a_n[inn]", 1);
  714.  
  715.   }
  716.  
  717.   if ((a_nbar = (double *) calloc (nnn_max + 1, sizeof (double))) == NULL)
  718.       fail_alloc ("a_nbar", 1);
  719.  
  720.   if ((sig_an = (double *) calloc (nnn_max + 1, sizeof (double))) == NULL)
  721.       fail_alloc ("sig_an", 1);
  722.  
  723.   if ((ndn = (int *) calloc (nnn_max + 1, sizeof (int))) == NULL)
  724.       fail_alloc ("ndn", 1);
  725.  
  726.   ma_tdt = p_of_nA (ndn, a_n, ns, vsa);
  727.  
  728.   gprintf (fpOut, "\nmean bubble domain area (. tdt): %lf\n", ma_tdt);
  729.   gprintf (fpOut, "domains grouped by coord no n=2(1)%d:\n", nnn_max);
  730.   for (inn = 2; inn <= nnn_max; inn++) {
  731.     for (idn = 0; idn < ndn[inn]; idn++) {
  732.       if (SHOW_ALL == 1) {
  733.         gprintf (fpOut, "...%5u ", *(a_n[inn] + idn));
  734.         if ((idn + 1) % 8 == 0)
  735.           gprintf (fpOut, "\n");
  736.       }
  737.       *(a_nbar + inn) += (double) (*(a_n[inn] + idn));
  738.     }
  739.     if (idn > 1) {
  740.       *(a_nbar + inn) /= (double) idn;
  741.  
  742.       gprintf (fpOut, "\n nof NN: %2d -- <A_n> ", inn);
  743.       gprintf (fpOut, "(mean of %d): %6.1lf \n",
  744.                idn, *(a_nbar + inn));
  745.       gprintf (fpOut, "\t -- <A_n>/<A>: %lf\n",
  746.                (double) (*(a_nbar + inn)) / ma_tdt);
  747.  
  748.       /* standard deviation */
  749.       for (idn = 0; idn < ndn[inn]; idn++) {
  750.         del_an = *(a_n[inn] + idn) - *(a_nbar + inn);
  751.         *(sig_an + inn) += del_an * del_an;
  752.       }
  753.       *(sig_an + inn) = sqrt (*(sig_an + inn) / (double) (idn - 1));
  754.       gprintf (fpOut, "\t  -- sigma_an: %lf\n\n",
  755.                (double) (*(sig_an + inn)) / ma_tdt);
  756.     }
  757.   }
  758.  
  759.  
  760. /*
  761.  * log vPoly area vs bubble area
  762.  */
  763.   gprintf (fpOut, "\n...bubble domain area vs Voronoi polygon area\n");
  764.   if ((dpa = (struct Pix *) calloc (ns, sizeof (struct Pix))) == NULL)
  765.       fail_alloc ("dpa", 1);
  766.  
  767.   nas = a_vs_vPa (dpa, ns, vsa, &ma_tdt, &mva);
  768.   gprintf (fpOut, "\tmean Vor poly area: %6.1lf:\n", mva);
  769.   gprintf (fpOut, "\t  mean bubble area: %6.1lf:", ma_tdt);
  770.  
  771.   gprintf (fpOut, "\t(%d unflagged Sites only)\n", nas);
  772.   for (is = 0; is < ns; is++) {
  773.     if (((dpa + is)->x > -1.0) && ((dpa + is)->y > -1.0)) {
  774.       gprintf (fpOut, "...Site %3d: %f %f\n",
  775.                (vsa + is)->sitenbr, (dpa + is)->x, (dpa + is)->y);
  776.     }
  777.   }
  778.  
  779.  
  780. /*
  781.  * boundary separation between NN domains: s_ij = d_ij - r_i - r_j
  782.  * note: each s_ij will occur twice
  783.  */
  784.   if ((bs = (float **) calloc (ns, sizeof (float *))) == NULL)
  785.       fail_alloc ("bs", 1);
  786.  
  787.   for (is = 0; is < ns; is++) {
  788.     if ((bs[is] = (float *) calloc ((vsa + is)->nnn, sizeof (float))) == NULL)
  789.         fail_alloc ("bs[is]", 1);
  790.  
  791.   }
  792.  
  793.   nbs = eval_sij (bs, ns, vsa);
  794.   gprintf (fpOut, "...%d values of pairwise NN bdy separations\n", nbs);
  795.  
  796.   for (is = 0; is < ns; is++) {
  797.     printf ("site %3d: ", (vsa + is)->sitenbr);
  798.     for (inn = 0; inn < (vsa + is)->nnn; inn++) {
  799.       printf ("...%6.1f ", *(bs[is] + inn));
  800.       if ((inn + 1) % 10 == 0)
  801.         printf ("\n");
  802.     }
  803.     printf ("\n");
  804.   }
  805.  
  806. /*
  807.  * evaluate histogram of boundary separations, normalized to 
  808.  * the mean bubble radius, <R> = sqrt( <A>/PI );
  809.  * -->note: each s_ij is counted twice
  810.  */
  811.   nh = 0;
  812.   for (is = 0; is < ns; is++)
  813.     nh += (vsa + is)->nnn;
  814.  
  815.   if ((hbs = (float *) calloc (nh, sizeof (float *))) == NULL)
  816.       fail_alloc ("hbs", 1);
  817.  
  818.   if ((ha = (float *) calloc ((size_t) N_BINS, sizeof (float))) == NULL)
  819.       fail_alloc ("hist\n", 1);
  820.  
  821. /*
  822.  * copy and normalize bdy separation values, eliminating entries 
  823.  * of 0.0 indicating that corresponding ref Site is out of bounds
  824.  */
  825.  
  826.   ih = 0;
  827.   for (is = 0; is < ns; is++) {
  828.     for (inn = 0; inn < (vsa + is)->nnn; inn++) {
  829.       if (*(bs[is] + inn) > 0.0) {
  830.         *(hbs + ih) = *(bs[is] + inn);
  831.         ih++;
  832.       }
  833.     }
  834.   }
  835.   printf ("\n...nbs: %d, nh: %d, ih: %d\n", nbs, nh, ih);
  836.  
  837. /*
  838.  * realloc
  839.  */
  840.   if (ih < nh) {
  841.     if ((hbs = (float *) realloc (hbs, ih * sizeof (float *))) == NULL)
  842.         fail_alloc ("hbs (realloc)", 1);
  843.   }
  844.   nh = ih;
  845.  
  846. /*
  847.  * evaluate statistics
  848.  */
  849.   gprintf (fpOut, "\n...statistics for boundary separation data\n");
  850.   mean = find_mean (hbs, nh);
  851.   sdev = find_sigma (hbs, nh, mean);
  852.   mr_tdt = sqrt (ma_tdt / PI);
  853.   gprintf (fpOut, "\tmean: %lf, mean/<R>: %lf\n", mean, mean / mr_tdt);
  854.   gprintf (fpOut, "\tsdev: %lf, sdev/<R>: %lf\n", sdev, sdev / mr_tdt);
  855.   gprintf (fpOut, "\t(in units of mean radius, <R> = %lf)\n", mr_tdt);
  856.  
  857. /*
  858.  * prepare to evaluate histogram: normalize and sort
  859.  */
  860.   for (ih = 0; ih < nh; ih++)
  861.     *(hbs + ih) /= (float) mean;
  862. #if defined(LINUX)
  863.   qsort (hbs, nh, sizeof (float), (__compar_fn_t) compare);
  864. #else
  865.   qsort (hbs, nh, sizeof (float), compare);
  866. #endif
  867.  
  868. #ifdef SHOW_HSORT
  869.   printf ("\n...sorted data:\n");
  870.   for (ih = 0; ih < nh; ih++)
  871.     printf (" hbs[%d] = %f\n", ih, *(hbs + ih));
  872. #endif
  873.  
  874.   HIST_HEADER = BS_HIST;
  875.   if (MIN_SET == 0)
  876.     min = (float) NSIJ_MIN;
  877.   if (MAX_SET == 0)
  878.     max = (float) NSIJ_MAX;
  879.   if ((-1.0e-10 < (min - max)) && ((min - max) < 1.0e-10))
  880.     min = max;
  881.   bw = (max - min) / ((float) (N_BINS - 1));
  882.  
  883.   printf ("\n...construct histogram of %d entries...\n", nh);
  884.   if (bw != 0.0)
  885.     construct_shist (nh, hbs, ha, bw, min);
  886.  
  887.   gprintf (fpOut, "%s\n", HIST_HEADER);
  888.   gprintf (fpOut, "...min/mean: %f, bw: %f, max/mean: %f\n", min, bw, max);
  889.   gprintf (fpOut, "...fractions\n");
  890.   for (i = 0; i < N_BINS; i++) {
  891.     gprintf (fpOut, " %6.2f", *(ha + i) / (float) nbs);
  892. /*      gprintf(fpOut, " %6.2f", *(ha+i) ); */
  893.     if ((i + 1) % 10 == 0)
  894.       gprintf (fpOut, "\n");
  895.   }
  896.  
  897.  
  898.   free (a);
  899.   free (area);
  900.   free (xs);
  901.   free (ys);
  902.   free (sig);
  903.  
  904.   for (inn = 0; inn < nnn_max; inn++)
  905.     free (a_n[inn]);
  906.   free (a_n);
  907.   free (a_nbar);
  908.   free (sig_an);
  909.  
  910.   free (dpa);
  911.  
  912.   for (is = 0; is < ns; is++)
  913.     free (bs[is]);
  914.   free (bs);
  915.   free (hbs);
  916.   free (ha);
  917.  
  918.  
  919. /*
  920.  * spatial correlations in the field A = A(r):
  921.  */
  922.  
  923. /*
  924.  * evaluate the average number, m(n), of NN of the NN of an n-fold
  925.  * coordinated site (->Aboav-Weaire law for cellular structures)
  926.  */
  927.  
  928.   if ((m_n = (float **) calloc (nnn_max + 1, sizeof (float *))) == NULL)
  929.       fail_alloc ("m_n", 1);
  930.  
  931.   for (inn = 0; inn < nnn_max; inn++) {
  932.     if ((m_n[inn] = (float *) calloc (ondn[inn] + 1, sizeof (float))) == NULL)
  933.         fail_alloc ("m_n[inn]", 1);
  934.  
  935.   }
  936.   if ((m_nbar = (double *) calloc (nnn_max + 1, sizeof (double))) == NULL)
  937.       fail_alloc ("m_nbar", 1);
  938.  
  939.   if ((sig_mn = (double *) calloc (nnn_max + 1, sizeof (double))) == NULL)
  940.       fail_alloc ("sig_mn", 1);
  941.  
  942.   for (inn = 0; inn <= nnn_max; inn++)
  943.     *(ndn + inn) = 0;
  944.  
  945.   nas = m_of_n (ndn, m_n, ns, vsa);
  946.  
  947.   gprintf (fpOut, "\n...average nof NN of the NN of n-fold sites\n");
  948.   gprintf (fpOut, "\t grouped by coord no n=2(1)%d:\n", nnn_max);
  949.   for (inn = 2; inn <= nnn_max; inn++) {
  950.     for (idn = 0; idn < ndn[inn]; idn++) {
  951.       if (SHOW_ALL == 1) {
  952.         gprintf (fpOut, "...%5.2f ", *(m_n[inn] + idn));
  953.         if ((idn + 1) % 8 == 0)
  954.           gprintf (fpOut, "\n");
  955.       }
  956.       *(m_nbar + inn) += (double) (*(m_n[inn] + idn));
  957.     }
  958.     if (idn > 1) {
  959.       *(m_nbar + inn) /= (double) idn;
  960.  
  961.       gprintf (fpOut, "\n nof NN: %2d -- <m_n> ", inn);
  962.       gprintf (fpOut, "(mean of %d): %5.2lf ",
  963.                idn, *(m_nbar + inn));
  964.  
  965.       /* standard deviation */
  966.       for (idn = 0; idn < ndn[inn]; idn++) {
  967.         del_mn = *(m_n[inn] + idn) - *(m_nbar + inn);
  968.         *(sig_mn + inn) += del_mn * del_mn;
  969.       }
  970.       *(sig_mn + inn) = sqrt (*(sig_mn + inn) / (double) (idn - 1));
  971.       gprintf (fpOut, " -- sigma_mn: %lf\n\n",
  972.                (double) (*(sig_mn + inn)));
  973.     }
  974.   }
  975.  
  976.  
  977. /*
  978.  * given the area, A0, of a bubble domain at location (x0, y0),
  979.  * evaluate the average, <Ak>, of areas of bubbles in the k-th NN
  980.  * shell of bubble 0, as a function of k
  981.  */
  982.  
  983. /* 
  984.  * construct kNN shells (topology: see also: xvora)
  985.  */
  986.   if (EVAL_KNN == 1) {
  987.     if ((sk = (struct kNNShell *) calloc (ns, sizeof (struct kNNShell))) == NULL)
  988.         fail_alloc ("sk", 1);
  989.  
  990.  
  991.     for (is = 0; is < ns; is++) {
  992.       if (((sk + is)->aka = (double *) calloc (KNN_MAX, sizeof (double))) ==
  993.           NULL)
  994.           fail_alloc ("sk->aka", 1);
  995.     }
  996.     for (is = 0; is < ns; is++) {
  997.       if (((sk + is)->tctc = (double *) calloc (KNN_MAX, sizeof (double))) ==
  998.           NULL)
  999.           fail_alloc ("sk->tctc", 1);
  1000.     }
  1001.  
  1002.  
  1003.     if ((sha = (struct linklist *) calloc (KNN_MAX, sizeof (struct linklist)))
  1004.         == NULL)
  1005.         fail_alloc ("sha", 1);
  1006.  
  1007.     gprintf (fpOut, "...constr kNN shells up to k=%d\n", KNN_MAX);
  1008.     frms = construct_kNNs (ns, vsa, sha, sk);
  1009.     gprintf (fpOut, "...avrg fractal meas of set: %lf\n", frms);
  1010.  
  1011.     gprintf (fpOut, "\n...<A_k>/<A_(k-1)> (<A> = %lf):\n", ma_tdt);
  1012.     gprintf (fpOut, "\t(sorted in order of area of root Site)\n");
  1013.     gprintf (fpOut, "\t(unflagged Sites only)\n");
  1014.  
  1015. /*
  1016.  * sort array sk in order of increasing areas associated with the root Site
  1017.  *
  1018.  * note: sorting scrambles the original Site indices so that flagged
  1019.  *   Sites are no longer correctly identified;
  1020.  */
  1021. #if defined(LINUX)
  1022.     qsort (sk, (size_t) ns, sizeof (struct kNNShell), (__compar_fn_t) cmp_root_a);
  1023. #else
  1024.     qsort (sk, (size_t) ns, sizeof (struct kNNShell), cmp_root_a);
  1025. #endif
  1026.  
  1027.     for (is = 0; is < ns; is++) {
  1028.       gprintf (fpOut, "...Site %3d ", (sk + is)->sitenbr);
  1029.       gprintf (fpOut, "A: %6.2lf A/<A>: %6.2lf ",
  1030.                (sk + is)->aka[0], (sk + is)->aka[0] / ma_tdt);
  1031.  
  1032.       /* normalize <A>_kNN by average over prev shell */
  1033.       for (k = 1; k < (sk + is)->kNN; k++) {
  1034.         gprintf (fpOut, "...%6.2lf ",
  1035.                  (sk + is)->aka[k] / (sk + is)->aka[k - 1]);
  1036.  
  1037.         if ((k + 1) % 8 == 0)
  1038.           gprintf (fpOut, "\n");
  1039.       }
  1040.       gprintf (fpOut, "\n");
  1041.     }
  1042.  
  1043. #ifdef FRMS_DEBUG
  1044.     KNN_MAX = 0;
  1045.     for (is = 0; is < ns; is++) {
  1046.       printf ("Site %3d: ", (vsa + is)->sitenbr);
  1047.       if (((vsa + is)->eFlag == UnSet) &&
  1048.           ((vsa + is)->aoiFlag == UnSet)) {
  1049.         printf ("kNN: %3d, frms: %6.3lf\n",
  1050.                 (sk + is)->kNN, (sk + is)->frms);
  1051.  
  1052.         if ((sk + is)->kNN > KNN_MAX)
  1053.           KNN_MAX = (sk + is)->kNN;
  1054.       }
  1055.       else {                    /* Site flagged out-of-bounds */
  1056.         if ((sk + is)->kNN == 0)
  1057.           printf (" out of bounds: -->ignored\n");
  1058.         else
  1059.           printf (" something wrong\n");
  1060.       }
  1061.       printf ("\tactual KNN_MAX: %3d\n", KNN_MAX);
  1062.     }
  1063. #endif
  1064.  
  1065.  
  1066. #ifdef KNNS_ARRAY_VER
  1067.     for (is = 0; is < ns; is++) {
  1068.       printf ("Site %3d\n", (vsa + is)->sitenbr);
  1069.       for (k = 0; k < (sk + is)->kNN; k++)
  1070.         show_kNN_list (&((sk + is)->sha[k]), k);
  1071.     }
  1072. #endif
  1073.   }
  1074.  
  1075.  
  1076. /* dealloc resources */
  1077.   fclose (fpOut);
  1078.  
  1079.   free (vsa);
  1080.  
  1081.   for (inn = 0; inn <= nnn_max; inn++)
  1082.     free (m_n[inn]);
  1083.   free (m_n);
  1084.   free (m_nbar);
  1085.   free (sig_mn);
  1086.  
  1087.   if (EVAL_KNN == 1) {
  1088.     for (is = 0; is <= ns; is++)
  1089.       free ((sk + is)->aka);
  1090.     for (is = 0; is <= ns; is++)
  1091.       free ((sk + is)->tctc);
  1092.     free (sk);
  1093.     free (sha);
  1094.   }
  1095.  
  1096.   free (ondn);
  1097.   free (ndn);
  1098.   if (CONTRACT == 1)
  1099.     free (cvsa);
  1100. }
  1101.